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I. INTRODUCTION 

Nonlinear dynamics of Josephson junction arrays 
(JJAs) has been a subject of extensive experimental and 
theoretical research [TJ [2]. The resistively and capaci- 
tively shunted junction (RCSJ) model of these arrays is 
described by the discrete sine-Gordon (DSG) equation 
which is ubiquitous in nonlinear physics [3l |4] . Among 
the actively discussed problems for the JJA dynamics the 
problem of the topological soliton (fluxon) response to 
the time-periodic bias, including the fluxon depinning, 
remains to be important. Properties of the small rf- 
biased Josephson junctions have been extensively studied 
both experimentally (starting from the pioneering papers 
of Shapiro [S]) and theoretically (with the focus on the 
phase- locking [6] and chaotic regimes [3|8]). In partic- 
ular, the rf-biased Josephson junctions have been used 
as a voltage standard (SKQ] . It is well-known [51 13] that 
a fluxon in a JJA is pinned due to the Peierls-Nabarro 
(FN) potential, unless a sufficiently strong bias is ap- 
plied. While the depinning of nonlinear excitations un- 
der the dc bias has been studied relatively well (see, e.g., 
Refs. [21 [T0Hl2] ). the problem of the ac depinning re- 
mains scarcely studied. 

The problem of ac fluxon depinning can be split into 
two cases: the symmetric and asymmetric depinning. 
The former case normally corresponds to the single har- 
monic driven systems, particularly being investigated in 
connection to the domain wall depinning in disordered 
systems [T31[I1]. In the latter case, a temporarily asym- 
metric (but with zero mean value) ac drive is applied. 
Here asymmetric drive means that the symmetry of the 
driving function is lowered, for instance, by applying 
a biharmonic signal. More details will be given in the 
next section. This phenomenon is based on the so-called 
ratchet effect, which is a unidirectional unbiased trans- 
port induced by symmetry breaking and nonlinearity |15l - 
117] . The mechanism of the phenomenon is based on the 
breaking of the symmetries connecting orbits with op- 
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posite velocities in the phase space [THl [19] and on the 
phase locking of the particle dynamics to the external 
drive [5D1[5T]. The rectification due to unbiased and tem- 
porarily asymmetric drive (normally biharmonic, con- 
sisting of a sinusoidal signal and its overtone) has been 
studied theoretically |27 and experimentally 23 in rf- 
biased small Josephson junctions. Also, the biharmonic 
drive has been used for chaos supression (see, e.g., Refs. 
[24l [25]). The experimental observation of the fluxon 
ratchet in Josephson junctions has been reported in sev- 
eral papers. Thus, fluxon ratchet in long Josephson junc- 
tions (LJJ) embedded in the inhomogeneous magnetic 
field that emulates spatially asymmetric sawtooth ratchet 
potential have been observed [IB]. In Ref. [57] the spa- 
tial asymmetry has been achieved by injecting external 
bias in a way that an effective asymmetric potential is 
created. The fluxon ratchet has also been observed in 
a JJA where spatial symmetry was broken by introduc- 
ing spatial modulation of the coupling energy [3H]. An 
interesting application of the fluxon ratchet effect as a 
fluxon pump has been suggested in Ref. (29j. Finally, 
the experimental observation of the fluxon ratchet in an 
annular LJJ due to temporal asymmetry (biharmonic rf 
bias) has been reported [S^. It should be emphasized, 
that the fluxon ratchet due to temporal asymmetry has 
not yet been investigated experimentally in discrete sys- 
tems. 

The ratchet phenomena have already been studied the- 
oretically for the nonlinear excitations in different dis- 
crete media including JJAs for both the spatial [5T| - f55] 
and temporal symmetry |36H38| breakings. In particu- 
lar, the difference of the ratchet motion in discrete and 
continuous systems has been highlighted (e.g., the dev- 
ils staircases and the non-zero depinning threshold 36 ). 
Nevertheless, the transition from a pinned to a running 
state in ac-driven ratchets requires a special investiga- 
tion. Also, it is worthwhile to note that finding mobile 
fiuxons in the limit of small coupling constants (k <^ 1) is 
an interesting challenge on its own. Therefore, address- 
ing the above questions is the main aim of this paper. 

The paper is organized as follows. In Section [llj we 
present the model. In the next section it is explained 
how the mode-locked solutions are computed and how 
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their properties are investigated. Linear stabihty of the 
single-harmonicany driven J JA is briefly discussed in Sec- 
tion |IV[ Asymmetric fluxon depinning under the bihar- 
monic drive is described in Section|V] FinaUy, Section [VT| 
contains the main conclusions of the paper. 



II. THE MODEL 

In this paper, we discuss the dynamics of an ar- 
ray of parallelly shunted and ac-biased small Joseph- 
son junctions within the RCSJ model. This system is 
described by the ac-driven and damped discrete sine- 
Gordon (DSG) equation, which can be written in a di- 
mensionless form as follows 

(pn — K A(/)„ + sin (I)n + a4>n + E{t) — 0, n — 1,2, . . . , N . 

(1) 

Here (jjn corresponds to the phase difference of the wave 
functions at the nth junction, A0„ = (/)„+i — 2(/)„ -I- 0„_i 
is the discrete Laplacian. The coupling constant k — 
•y/$o/(27r/cL) measures the discreteness of the array, 
where $o is the magnetic flux quantum, L is induc- 
tance of an elementary cell, and Ic is the critical cur- 
rent of an individual junction. The dimensionless dissi- 
pation parameter is then a = $0/ (27r Jci?) , where R is 
the resistance of an individual junction, and the time is 
normalized to the inverse Josephson plasma frequency 
l/ujQ = \/ C'^i^)/ {2i:Ic) with C being the junction capaci- 
tance. Finally, E{t) — E{t + T), T = 2'k/uj is an external 
bias, which has a zero mean value [{E{t))t = 0], applied 
to each junction of the array. In the following, we assume 
E{t) to be of the form 

E{t) ^ EiCOs{ujt) + E2Cos{2u!t + e), (2) 

Notice that the superposition of the two harmonics makes 
the periodic force to be asymmetric in time for almost 
all values of 0, a feature which can be used to break the 
temporal symmetry of the system (more details will be 
given below). Only the circular arrays with are to be 
considered, therefore the periodic boundary conditions 
apply: 0„+Ar(i) = (/)„(<) -l- 2Qtt, (pn+Nit) = <j)„{t), where 
Q is an integer constant that stands for the net number 
of fluxons trapped in the array. Further on only the case 
of one fluxon {Q = 1) will be considered. 

In JJAs the topological solitary waves have the phys- 
ical meaning of trapped magnetic flux quanta (fluxons) 
and the average voltage drop reads 

. N t 
^=^E 1™ 7 / ■ (3) 

iV ^ — ' t-i-oo t /n 

If the fluxon is moving with a non-zero net velocity then 
V ^ and V = otherwise. 

The experiments with annular JJAs have been per- 
formed for typical lengths N - 8-=-30 (see Refs. [Ttl^l^). 
In the following, we consider the case of an array with 



N = 10 junctions subjected to periodic boundary condi- 
tions (annular array) unless stated otherwise. 

The unidirectional fluxon motion can take place either 
on regular trajectories (limit cycles) or on chaotic trajec- 
tories. Further on we will refer to the trajectories with 

7^ as to the transporting ones, while the trajectories 
with = will called non-transporting trajectories. Ob- 
viously, only the transporting trajectories are of interest 
within this paper. 

The regular transporting trajectories correspond to the 
limit cycles of Eq. ([!]) , which are mode-locked to the fre- 
quency of the external bias. On this orbit, the average 
kink velocity is expressed as (v) = j • 5^ , where the wind- 
ing numbers k and I are integer. In the resonant regime, 
the fluxon travels k sites during the time IT = 2ttI/uj, 
so that, except for a shift in space, its proflle is com- 
pletely reproduced after this time interval (in the pen- 
dulum analogy, this orbit corresponds to k full rotations 
of the pendulum during I periods of the external drive). 
The voltage drop in an annular JJA (with one fluxon in 
it) is related to the average fluxon velocity {v) by the 
equation V = 2tt{v)/N [5]. 

According to Ref. [33] i in order to obtain the directed 
fluxon motion, all symmetries that relate two fluxons 
with opposite velocities should be broken. It happens if 
the following inequality is true: E{t) ^ -E{t -I- T/2). 
The bias (pi) satisfles it, but in the strictly Hamilto- 
nian case (a — 0) an additional condition is necessary: 
E{-t + t') ^ E{t). 

III. FLOQUET ANALYSIS OF THE 
MODE-LOCKED STATES 

In order to understand better the depinning transition, 
we focus first on the regular mode-locked solutions (limit 
cycles) of the driven DSG equation. The fluxon periodic 
orbit is computed by flnding zeroes of the map 

ife,(r)x = x, (4) 

where the vector X consists of the dynamical variables 
{0n, <i}n}n=i- The Operator Iki stands for the integration 
of the equations of motion ([T]) during the time IT and 
afterwards the shift of the final solution by k sites forward 
if fc < or backward if fc > 0. The case /c = corresponds 
to the fluxon pinned to a lattice site. 

A flxed point of the map Q is a mode-locked solu- 
tion {(j)n\t), (l)n\t)}n=i which reproduces itself after the 
time IT with the space shift by k lattice sites backward 
or forward. Next, we substitute the expansion 

(/.„(t) = (/.i°)(<)+£„(t) , (5) 

into Eq. ([T]). For the case of stondm^ fluxon (fc — 0) after 
keeping only the linear terms, we obtain the following set 
of linear ODEs with periodic coefficients: 

e„ = -Q;£„-l-KAe„-cos[(/)^°' {t)]e„ , n = 1,2, . . . ,N. (6) 
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The map 
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is constructed from the solutions of the system It 
relates the small perturbations e{t) = {£n{t)}n=i at the 
time moments i = and t = IT. The 2N x 2N Flo- 
quet (monodromy) matrix M contains all the necessary 
information about the linear stability of the system. If 
this matrix has at least one eigenvalue with |A„| > 1 
(n = 1,2, . . . ,2N), then the system is unstable. If for 
all eigenvalues |A„| < 1, the system is stable. It is well- 
known [3D] that these eigenvalues come in quadruples, so 
that if A„ is an eigenvalue, then A*, i?/A„ and i?/A* 
(here R = e"'"'"/", see, for example, Refs. gllSl]) are 
also eigenvalues. Thus, the Floquet multipliers lie either 
on the circle of the radius R (will be referred to as a R- 
circle throughout the paper) or may depart from it after 
collisions. The notable difference of the ac-driven case 
from the dc-driven (autonomous) case is the absence of 
the degeneration with respect to time shifts, which man- 
ifests itself in the absence of the eigenvalue A = 1 |43] . 

With the help of the Newton-Raphson iterative 
method it is possible to compute numerically the respec- 
tive mode-locked limit cycle for the given period IT with 
a desired computer precision. For details one might con- 
sult Ref. [44|. The advantage of this approach is that 
not only attractors, but also repellers, can be computed. 
Also, wrong conclusions which can be made due to sen- 
sitivity to initial conditions can be avoided. 

In this paper, we plan to compute the mode-locked 
limit cycle that corresponds to the standing fluxon and 
to path-follow it while a control parameter is changed un- 
til the cycle becomes unstable or completely disappears. 
By monitoring the moduli of the Floquet eigenvalues |A„| 
one can obtain the information about the underlying bi- 
furcations and, consequently, about the depinning tran- 
sition. 



IV. FLUXON DYNAMICS UNDER THE 
SINGLE-HARMONIC DRIVE 

Before investigating the problem of the directed fluxon 
motion under the influence of a biharmonic signal ([2| we 
briefly consider the case of a single harmonic drive when 
E2 — 0. It is of interest to investigate the stability of a 
mode-locked state that corresponds to a standing fluxon. 
The existence of the non-zero Peierls-Nabarro (PN) po- 
tential causes fluxon pinning to the lattice [45]. Intu- 
itively, it is not hard to understand that on the parameter 
plane (k, Ei) one can draw a curve that separates the area 
where only the stable mode-locked standing fluxons ex- 
ist (the pinning area) from the are where mobile fluxons, 
both chaotic and regular, coexist together with the stand- 
ing fluxons (the transporting area). This curve should 
exist as a dependence e['^'^ = E^^'^ (k) where for a given 



K only pinned mode-locked states exist if Ei < E' 

El > e[''\ the dynamics appears to be more complex 
including diffusively moving fluxons, mode-locked mov- 
ing fluxons or chaotic dynamics of the whole array when 
individual fluxons cannot be identified. The latter case 
corresponds to the non-existing area and will not be dis- 
cussed in this paper. In the limit k — >■ 00, the effects 

(c) 

of discreteness disappear, thus E^ should decrease. On 
the other hand, the decrease of k means that the PN 
barrier becomes stronger and thus a larger amplitude is 

(c) 

necessary to overcome it. As a result, E^ increases when 
K — ^ 0. The exit from the pinning area [below the curve 

isj^-* (k)] can lead to different scenarios depending on the 
direction of the exit. 

The issue of discrete kink unlocking (depinning) in the 
DSG lattice has been studied in Ref. [35] j but for the 
parameter ranges which are different from ours. We are 
dealing with large bias amplitudes and small couplings. 
Especially it is useful to focus on the Floquet analysis 
of pinned fluxon states in the limit of small (k ^ 1) 
couplings. Also, this section will make more clear the 
presentation of the original results in the next sections. 

The behavior of the Floquet eigenvalues (computed as 
described in Sec. Ill) is shown in detail in Fig. [I] In 
panels (a)-(b), their evolution as a function of the cou- 
pling constant k is given starting from the anticontinuum 
(k — 0) limit. Note that the eigenvalues are placed on 
the complex plane symmetrically with respect to the real 
axis, thus we limit ourselves with the eigenvalues lying in 
the upper half-plane (Re A > 0). The bias amplitude now 
is El = 0.05, so that the excitations on the fluxon back- 
ground can be treated as small. In the anticontinuum 
limit, all the eigenvalues sit in one point and with the 
growth of K they separate forming two distinct groups: 
the modes associated with the linear spectrum (Joseph- 
son plasmons) and the internal mode(s). The plasmon 
band extends linearly with k according to the dispersion 
law 



wl(9) = -1/1 -I- 4k sin 



2 Q 



(8) 



Due to finiteness of the array, the wavenumber q G [0, 2tt) 
attains only discrete set of values qm = 2'Km/N , m = 
±.\,...,±N. The localized eigenmode detaches itself 
from the linear band: as it is well known in the discrete 
Klein-Gordon theory [37], the internal mode is soft, so 
that it detaches from the linear band into the gap. In 
panels (c)-(l) of Fig. [l| the shapes of the eigenvectors e 
(both real and imaginary parts) are shown for the param- 
eter values El — 0.05 and n = 0.1. These eigenvectors are 
ordered from top to bottom according to the decrease of 
argA„. From the shape of the eigenvectors one can con- 
clude that all of them, except one, are delocalized and 
therefore they are associated with the linear spectrum. 
The eigenmode in Fig. [ijf) has a clear localized struc- 
ture. In our case, the fluxon state is locked to the external 
drive with the frequency uj = 0.35 which lies in the gap 



4 



arg A 




(e) 

(f) 
(g) 



argA 



(k) 



(1) 




|A| 



G 0.57S 0.5S0 0.5S2 





FIG. 1: (Color online) Phases (a) and moduli (b) of Floquet eigenvalues as function of k for Ei = 0.05, a = 0.1, and lj — 0.35. 
Floquet eigenvectors [panels (c)-(l)] at Ei = 0.05, k = 0.1 (see text for details). Panels (m)-(p) show the phases and moduli 
of Floquet eigenvalues at k = 0.1 as function of E\. Upper inset in panel (b) shows the eigenvector for unstable eigenvalue at 
K = 0.5058245 [Re £„ (+) and Im e„ (♦)]. Same for the inset in panel (n) at E\ — 0.58183325. The lower insets shows the 
positions of eigenvalues. 



of the spectrum ([8]). Only the overtones n > 3 can lie in 
the linear spectrum: 1 < nui < vT^Mv?, therefore the 
excited linear modes satisfy the condition nu — ujL{Qm), 
TO = l,2,...,iV. The increase of k widens the linear spec- 
trum and, as a result, the respective eigenvalues spread 
around the i?-circle. The collisions of the eigenvalues at 
argA — TT with the highest (m = N) linear mode are 
marked by the vertical red lines at k = [(7aj/2)^ — l]/4 
and K — [(9a;/2)^ — l]/4. They correspond to the para- 
metric resonances of this mode with the external drive: 
^l{(1n) = 7a;/2 and ujl{<1n) = 9i^/2. Collisions at 
argA = marked by the blue lines at k = [(Sw)^ — l]/4 
and K = [(4ci;)^ — l]/4 and correspond to the main res- 
onances Wi(q'jv) — 4w and u)L{qiq) ~ Alj. The internal 
mode evolution with varying k, can be clearly seen due to 
its nonlinear behavior as a function of k, as compared to 
the linear evolution of plasmon modes. After interactions 
with the linear spectrum (marked by the avoided cross- 
ings), the internal mode hits the real axis at argA = 0, 
thus signaling the resonance with some overtone of the 
driving frequency. After that the eigenvalue increases ex- 
tremely fast [note the almost vertical growth of |A„| in 
Fig. [ijb)], leaves the i?-circle and finally exits the unit 
circle. We conclude that the disappearance of the stand- 
ing mode-locked fluxon takes place via the tangential 
(saddle-node) bifurcation. The excitation of the fluxon 
along the destabilizing direction in the phase space leads 



to chaotic fluxon diffusion with the net zero velocity. 

While the previous investigation took place at rather 
small amplitudes for which the excitation on the fluxon 
background can be considered as small, now we focus on 
the evolution of the Floquet eigenvalues as a function 
of the driving amplitude for the fixed value of coupling 
K = 0.1. As one can see in Fig. [T]jm-n), the eigen- 
value that corresponds to the localized mode moves on 
the i?-circle in the direction opposite to the motion of 
the linear band. After following the respective eigen- 
value as a function of -Ei, we see that it collides with 
its complex conjugated counterpart at argA = tt (see 
the resulting "bubble" in the |A(i<^i)| dependence), con- 
tinues its motion along the JJ-circle, then interacts with 
the eigenmodes of the linear spectrum and finally collides 
with its complex conjugate counterpart on the real axis 
but at argA — 0. Next, the modulus of this eigenvalue 
grows fast and finally exceeds the unit circle. Again, 
the standing mode-locked fluxon state disappears via the 
tangential bifurcation. The unstable eigenvector is again 
spatially localized (see the inset in Fig. [l]i). The pertur- 
bation of the kink in the unstable direction leads to the 
growth of a localized oscillation on the kink background 
that results in the kink-antikink pair birth and eventual 
destruction of the one-soliton state. Other scenarios are 
also possible. Similar investigation for the case k, = 0.2 
shows that the system can undergo additional bifurca- 
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tions and a new standing mode locked fluxon can appear 
with a slightly different shape. However, the final disap- 
pearance of the standing kink goes similarly to the above 
scenario via the tangential bifurcation. In this case, the 
time evolution of the unstable fluxon turns into chaotic 
diffusive jumps of the kink backward and forward with 
the net zero velocity. 

The main conclusion of this subsection is that stand- 
ing mode-locked fluxon states disappear via a tangential 
bifurcation, when the driving amplitude Ei or the cou- 
pling constant exceeds some critical value. The resulting 
dynamics can either lead to fluxon chaotic diffusion or to 
chaotic motion of the whole lattice and consequent fluxon 
destruction. Our analysis with TV = 20 and iV = 30 has 
shown that the size effects are insignificant: the critical 
values at which the tangential bifurcation takes place dif- 
fer very slightly. The instability evolves via the excitation 
of a localized internal mode. This is in accordance with 
the previous results [36, 42^ on the internal mode role 
in the kink mobility. More details on the kink response 
to the single harmonic ac drive (including the kink dy- 
namics) can be found in Ref. |46| . In particular, it has 
been shown that in the unlocking limit the kink dynamics 
shows the type-I interniittency. 



V. ASYMMETRIC FLUXON DEPINNING 

Now we focus on the main task of this paper, namely 
on the studies of the fiuxon motion when the bias ^ 
is biharmonic {E2 7^ 0). Here and further on, we take 
El = E2 unless stated otherwise. According to the pre- 
vious studies |331 136] , the fluxon motion becomes unidi- 
rectional due to the symmetry breaking and its average 
velocity depends solely on the system parameters. Also, 
from those papers we already know that critical parame- 
ter values should exist above which the applied external 
bias overcomes the forces of pinning and the fluxon starts 
to propagate. In particular, such critical values exist for 
the amplitudes of the external drive E12 and the phase 
shift 0. 

Including the second harmonic in the bias ([2|, we add 
the third parameter (the phase shift 9) into consider- 
ation and consequently we must think about the exis- 
tence diagram in the three-dimensional parameter space 
(k, El, 9). However, since the average voltage drop satis- 
fles V{9) — V{9 + 2tt), it is possible to construct the exis- 
tence diagram on the plane (k, Ei). This diagram for the 
existence of moving fluxon states is shown in Fig. [2j The 
numerical experiment has been performed as follows: the 
pinned state with the fluxon oscillating around its center 
of mass has been continued with small increase of Ei or 
K until this state loses stability and starts to propagate. 
Transporting trajectories are mostly chaotic, but regular 
mode-locked trajectories may also exist. These trajec- 
tories will be analyzed in detail in the next subsections. 
As one can see, similarly to the case of the single har- 
monic drive, there exists a curve e[''^ = E^'^ (k) ( shown 



by markers) that separates the pinned and transporting 
parts of the diagram. In the area down and to the left 
from the curve i?^^-' (k), there are only mode- locked stand- 
ing fluxons and no mobile fluxons exist {V = 0) for any 
9 e [0, 27r). The area to the right side of the curve cor- 
responds to the situation where moving fluxons {V 7^ 0) 




K 

FIG. 2: (Color online). Existence diagram of moving fluxons 
on the plane (k, £1) at a = 0.1 for u = 0.25 (e), cj = 0.35 (♦) 
and ui = 0.5 ( x ) . Solid lines are used as a guide for an eye (see 
text for details). The upper and lower insets show the time 
evolution of the central junction for 6 = 2.5, uj — 0.35, a = 0.1 
and K = 0.12, Ei = E2 ^ 0.22 (lower inset) and k = 0.09, 
El = E2 = 0.27 (upper inset). Different colors correspond to 
different initial times: to = (black) and to = n/uj (blue). 

exist for some values of 9. Obviously, in the continuum 
limit K — >■ 00, the critical depinning amplitude tends to 
zero, while in the anticontinuum limit k — > 0, it is neces- 
sary to increase the driving amplitude in order to unlock 
the fluxon. At some point the critical drive is too strong 
and the fluxon is destructed due to the fluxon-antifluxon 
pair creation and the chaotic dynamics of the whole lat- 
tice. An example of such a pinned mode-locked trajec- 
tory is shown in the lower inset of Fig. [2j Note that the 
pinning to the lattice is independent on the initial con- 
ditions. For example, in this figure two trajectories with 
the initial times to = and to = T/2 are shown. Also, 
the scenario is the same for other values of initial time 
to e [0, T), the initial fluxon position, and different initial 
kicks. The examples of transporting chaotic trajectories 
with different initial times (<o = and to = 2^/2) are 
illustrated by the upper inset of Fig. [2] Independently 
of the initial conditions at t — >■ 00, the system settles 
on the chaotic attractor that corresponds to the directed 
fluxon motion with the same average velocity (in this case 
F~ 0.0055). 

Note that the principal characteristics of the e[''\k) 
dependence remain qualitatively the same for different 
values of Lo. It has been shown in Ref. [3^ that the 
ratchet fluxon motion is well pronounced in the fre- 
quency range < w < a;L(0)/2. The Peierls-Nabarro 
frequency (which equals the frequency of the internal 
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mode 'SO') at k = 0.1 (according to the relation ojpN — 
A/27raoK3/2 exp (-7r2^/2), ao ~ SOtt , see Ref. [5T] ) 
equals wp^r ~ 0.91. Thus, for small couplings the work- 
ing frequency range lies below all characteristic system 
frequencies. Among the other system parameters the 
coupling constant and dissipation depend on the phys- 
ical characteristics of the specific array and therefore 
they cannot be changed easily in the experiment. Con- 
sequently, we are left with the bias parameters Ei and 
which can be tuned at will. Therefore it would be of 
interest to study the depinning process as a function of 9 
(with other parameters fixed) and the same for Ei. The 
following subsections will be devoted to this task. 



A. Depinning as a function of the bias amplitude 

For continuous ratchets it was shown previously |52j 
that the average kink velocity is proportional to EfE2, 
so that, provided the respective symmetries are broken, 
the directed motion can occur for arbitrary small values 
of the driving amplitudes. For the DSG equation, the 
dependence of the average kink velocity on the driver 
amplitudes is not a smooth monotonic function of E12, 
but a complex function that either entirely or partially 
consists of plateaux of different lengths that correspond 
to mode-locked fluxon states. At k > 1 this dependence 
resembles closely a "devil's staircase", while this resem- 
blance disappears when k is decreased leaving only iso- 
lated islands of regular motion [36j . 

In Fig. [3] the dependence of the Floquet eigenvalues 
on the driving amplitude is shown for different values of 
K. Comparing Figs. [T|m-n) and[3]^a-b), we observe the 
close similarity in the Floquet eigenvalue behavior in the 
single harmonic and biharmonic biased cases. It appears 
also that independently on the coupling strength, the loss 
of stability of the limit cycle appears through a collision 
of two eigenvalues at argA — and the departure of 
one of the eigenvalues out of the unit circle [see the inset 
in Fig. [Sj^b)] that signals a tangential bifurcation tak- 
ing place at some critical value of the driving amplitude. 
This critical value decreases with the growth of k. Note 
that the critical value of the bias amplitude is approxi- 
mately two times smaller than in the single harmonic case 
because now the bias consists of two harmonics with the 
amplitude Ei. In panels (a-b), it should be noted that 
before the tangential bifurcation several period-doubling 
bifurcations take place that can easily be identified by 
the bubble-like deviations on the dependence |A„(i?i)| 
from the value |A„| = R. The amplitude of these devia- 
tions does not exceed the value |A| = 1 and they decrease 
with K, so that the limit cycle remains stable. At higher 
K = 0.2, the largest bubble corresponds to a pair of eigen- 
values colliding away from the real axis (that corresponds 
to the Hopf bifurcation). For even higher k = 0.5, again 
a bubble associated with the period-doubling bifurcation 
can be seen. Although for the given range of parameters 
these bifurcations do not cause instability of the limit 



cycle, the increase of w or decrease of a may change the 
nature of the depinning transition. The computations 
have been performed for different values of 9 and the de- 
pendencies A„(£^i) appear to be rather similar since the 
curves of different colors almost coincide. The only dif- 
ference lies in the different critical values of Ei , when the 
tangential bifurcation takes place, however these differ- 
ences are rather small. 

The unstable eigenvector, shown in the inset of panel 
(d) has a pronounced localized asymmetric shape. The 
eigenvectors that correspond to the opposite (with re- 
spect to the shift 9 — >■ 9 zt tt) values 9 = 1.5 and 
9 = 1.5 — TT ~ —1.64159 are related to each other ap- 
proximately by the inversion e — > — Ce, where C is an 
arbitrary constant that appears because Eqs. ^ are lin- 
ear. This means that the unstable perturbation drives 
the fluxon in the opposite directions for 9 and 9 + n, 
respectively. 

The behavior of the eigenvalues in the parameter space 
{k,Ei,9) has very complicated structure. In particular, 
it is possible that for some values of 6', the standing fluxon 
states can exist after the tangential bifurcation shown in 
Fig. [3] For the sake of brevity, we are going to refer 
to it as to the first tangential bifurcation. Indeed, such 
an example is given by Fig. |4j After the firts tangential 
bifurcation, the standing solution (limit cycle) becomes 
unstable at Ei ~ 0.2357 and soon disappears but reap- 
pears back after another tangential bifurcation. In the 
short interval between these bifurcations, a chaotic mov- 
ing fluxon appears. The further increase of the driving 
amplitude leads to the Hopf bifurcation (collision of two 
eigenvalues on the complex plane) and to the eventual 
departure of two eigenvalues out of the unit circle. This 
instability leads to formation of a large localized excita- 
tion on the fluxon and the eventual creation of fluxon- 
antifluxon pairs. Thus, it is possible to pass from the 
pinned area to the non-existing area (where no fiuxons 
exist due to chaos) omitting the transporting area (where 
the moving fiuxons exist). 

This bifurcation analysis explains the non-monotonous 
behavior of the average velocity dependence on the bias 
amplitude for small k (see Fig. 8 of Ref. [SB]), when 
the intervals of zero and non-zero fiuxons velocities in- 
terchange with the growth of Ei . 



B. Depinning as a function of the phase shift 

The directed fluxon transport in the continuous sine- 
Gordon model shows [55] that the dependence of the 
average fluxon velocity (and, consequently, the voltage 
drop) behaves as sm[9 — 6*0(0;, a)]. As it has been shown 
in Ref. in the weakly discrete case, the two values 
of 9 where V{9) = become intervals, and the size of 
these intervals increases when k decreases. Moreover, 
the dependence V{9) changes from the continuous be- 
havior into the complex piecewise function which may 
include resonant plateaux V = ujk/{lN) and loses its re- 
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FIG. 3: (Color online). Dependence of the phases (upper figures) and moduli (lower figures) of Floquet eigenvalues on bias 
amplitude _Ei at a = O.f, cj — 0.35, and for different coupling constants and 6: k — 0.09 (a-b), 8 — 1.5 (black lines), 6^ = 
(green lines); k = 0.2, (c-d), 6 = 1.5 (black lines), 61 = (green lines); 6 = -1.64159 (blue lines); k = 0.5, 6 = 1.5 (e-f). The 
inset in panel (b) shows the position of Floquet eigenvalues on complex plane at Ei = 0.04511641. The inset in panel (d) shows 
the profile of unstable eigenvectors Re£„ (+) and Ime„ (♦) at Ei = 0.1403205, 9 = 1.5 for eigenvalue A = 1.005497 (black 
lines) and Ei = 0.14032053, = -1.64159 (scaled by 2, red lines) for eigenvalue A = I.0005I5 [Ree„ (x) and Im£„ (o)]. Solid 
lines are used as guides for an eye. 



semblance with the sine function as long as k — ?> 0. For 
K > 1 there exist at least two intervals on the 9 axis 
which we denote as (0^,02 ) and (0^,0^) where V ^ 0. 
In particular, V{e) < if 6* e (0^,0-) and ¥{6) > if 
9 G {9^,92)- It appears, however, that for small k the 
dependence V{9) may take a more complicated shape (for 
example, see Fig. 8 of Ref. |36]). This also can clearly be 
seen in the main graph of Fig. [5] where the dependencies 
for V{9) for several values of k and for the fixed value of 
El — 0.27 are given. This figure in some sense is com- 
plementary to Fig. |2] because it explains in detail what 
is happening around the curve i?^'^' (k) when 9 changes. 
While at K = 0.09 only two transporting intervals where 
V =^ aie observed with 9^ ~ -0.69, 61^ ~ 0.23 and 
9^ ~ 2.47, 9^ ~ -2.92. For other values, namely, for 
K — 0.12 there can exist two additional transporting in- 
tervals. The transporting intervals (shown by vertical 
bars on the upper inset of Fig. [5]) decrease as k — >■ 
at El — const. The same behavior occurs if k is fixed 
and El decreased. It is worth to notice the existence 
of regular transporting limit cycles that correspond to 
the constant voltage plateaux V = ±0.035 for k = 0.12. 
Their existence has been obtained by the direct Runge- 



Kutta integration of the equations of motion and also 
has been confirmed by the Newton method with /c = ±1, 
1 = 1. The careful investigation of the transition between 
the pinned state and these plateaux shows the existence 
of narrow hystereses windows as shown in the lower in- 
set of Fig. [Sj In order to understand the nature of the 
additional pinning intervals, the Floquet analysis of the 
standing limit cycles has been performed as a function of 
9. 

In Fig. |6] the behavior of the Floquet eigenvalues as a 
function of the desymmetrizing parameter 9 for different 
values of Ei below and above the dcpinning threshold 
is shown. The left panels (a-b) correspond to the situa- 
tion when K = 0.09 and Ei is increased from below the 
critical depinning dependence E^'^\k). One can see that 
initially the curves arg A„(6') do not cross each other. For 
El = 0.26 all the eigenvalues lie on the i?-circle except for 
the two pairs with argA„ oscillating around ±37r/2 and 
having |A„| ^ R, as it can clearly be seen from Fig. |6][a- 
b). These eigenvalues correspond to the spatially local- 
ized eigenvectors as shown in Fig. 6f^c). Moreover, the 
non-localized eigenvalues are virtually independent on 9. 
The shape of the localized eigenvectors does not differ sig- 
nificantly for 9 = ±7r and 9 = as well. A slight increase 
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El El 



FIG. 4: Phases (a) and moduli (b) of Floquet eigenvalues as 
a function for uj — 0.35, a — 0.1, k = 0.12, 8 — 1.5. Panels 
(c) and (d) give details in interval Ei e (0.22,0.24). 
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FIG. 5: Average voltage drop as a function of phase shift 
6 for uj = 0.35, a = 0.1, Ei = E2 = 0.27, k = 0.09 (♦) and 
K = 0.12 (o). The upper inset shows the existence diagram of 
transporting trajectories (vertical bars) on plane {k,6). The 
lower inset shows the details of hysteresis (see text). 

of the driving amplitude to the value Ei — 0.265 leads to 
further increase of the oscillation of localized eigenvalues. 
Upon the further increase of Ei , the eigenvalue collisions 
at A = — 1 and out of the real axis take place (see, for 
example, the green line in Fig. [6]^ a) that corresponds to 
El = 0.269). The monitoring of the eigenvectors that 
correspond to the eigenvalues with |A„| 7^ R (see panel 
(d) of Fig. |6| shows that the eigenvector shape changes 
significantly, although they remain spatially localized. 
Now the phases of all the eigenvalues demonstrate sig- 
nificant dependence on 6. The crossings of the arg A„(0) 
curves demonstrate that the localized eigenmode "inter- 
acts" with the modes of the linear spectrum. First of all 
we note the collision of the eigenvalues at arg A = tt that 
corresponds to the period-doubling bifurcation. These 




FIG. 6: (Color online). Phases (a) and moduli (b) of eigen- 
values of Floquet matrix as a function of phase delay 9 
for UJ = 0.35, K = 0.09 and Ei = E2 = 0.26 (black). 
El = E2 = 0.265 (blue), Ei = E2 ^ 0.269 (green) and 
E\ — 0.27 (red). Panels (c)-(e) show the real and imagi- 
nary parts of destabilizing eigenvectors the eigenvalues that 
lie out of the ii-circle (see text for details): (c) Ei — 0.26 
with 6 = —TV [red, Re s (x), Im e (o)], 9 = — 7r/2 [blue. 
Re £ (*), Im e (o)] and e = [black,Re e (-I-), Im e (♦)]; (d) 
El = 0.269 (the same); (e) Ei = 0.27 with 9 = -2.914 [red. 
Re e (x), Im e (o)], 9 = — 7r/2 [blue. Re e (*), Im e (o)] and 
9 = -0.69206 [black, Re e (+), Im e (♦)]. 



bifurcations are clearly seen if one compares the lines for 
El = 0.265 and Ei = 0.269 in Fig. |6];a-b), however, they 
do not cause any instabilities. After reaching some criti- 
cal value of El , the localized eigenvalues which previously 
have resided in the complex plane collide with each other 
on the real axis at Re A > 0. The further increase of Ei 
leads to the fast exit of one of them out of the unit cir- 
cle and to the subsequent tangential bifurcation. Beyond 
this bifurcation point, the non-transporting limit cycle 
disappears. For Ei = 0.27 there exist four tangential bi- 
furcations, two direct and two inverse, that happen not 
far from the values 9 — and 9 — ±7r. The respective 
unstable eigenvectors are spatially localized, as shown in 
Fig.[6f^e). The respective repeller also has been computed 
(not shown in the figure for the sake of clearness), and 
its eigenvalues join the attractor eigenvalues at the bi- 
furcation point. The transporting intervals (0j~,0^) and 
{9i,92) in Fig. [5] for k = 0.09 obviously correspond to 
the windows in Fig. [6]ja-b) where the standing fluxon 
does not exist. Similarly to the dependencies shown in 
Fig. [3j the bifurcation points are preceded by the sharp 
singular-like growth of the modulus of the unstable eigen- 
value. 

A different scenario is presented in Fig. [7j As it has 
been discussed in the previous subsection, the existence 
subspaces of the transporting and non-transporting tra- 
jectories in the parameter space (k, Ei, 9) can take a pe- 
culiar and complicated shape. This is also true for the 
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discrete soliton ratchets with the spatial symmetry break- 
ing studied in Ref. ^33j (see Fig. 12 therein). If one stays 

above the critical dependence e[''\k), it is possible to 
observe the existence of a non-transporting limit cycle af- 
ter the first tangential bifurcation as shown in Fig. |4] It 
is logical to investigate such a case for different 6''s when 
El > e['^\k). This exactly has been done in Fig.jrjwhere 
the eigenvalue behavior has been plotted as a function of 
9 for El = 0.25 and Ei = 0.27. On these figures again 
four tangential bifurcations are seen (two direct and two 
inverse), clearly pointing out the windows where the di- 
rected fluxon motion takes place. These bifurcations are 
characterized by sharp singular growth of the respective 
|A„(6')| dependencies. 




FIG. 7: (Color online). Phases (a) and moduli (b) of eigen- 
values of Floquet matrix as a function of phase delay for 
K = 0.12, El ^ E2 = 0.25 (red) and Ei = E2 ^ 0.27. The 
rest of parameters is as in Fig. |6] The inset shows the unsta- 
ble eigenvector Re e {+) and Im e (♦) at 6 = —3.03 just after 
the period-doubling bifurcation (see text for details). 

Apart from these bifurcations which are typical for 
the case considered before (including Fig. |6]), the two 
large "bubbles" on the dependence |A„(6')| can be no- 
ticed. They correspond to the eigenvalue collisions on the 
real axis at Re A„ < and to the subsequent departure 
of one of them out of the unit circle. Thus, we observe 
period-doubling bifurcations that make the mode-locked 
state ustable. Along the interval tt < 9 < tt the four 
period-doubling bifurcations take place: two direct and 
two inverse. The destabilizing eigenvector is spatially lo- 
calized as in the case of tangential bifurcations (see the 
inset of Fig. [t]). Note that the "bubbles" of the insta- 
ble eigenvalues increase when Ei increases. The stable 
period-2 (with / = 2) limit cycles that appear after the bi- 



furcation, have been computed with the Newton method 
until the next period-doubling bifurcation takes place. 
The emerging period-4 cycle has been computed as well. 
In the next section, it will be shown that in fact the cas- 
cade of period-doubling bifurcations takes place leading 
to chaotic dynamics. 

Now we can explain why the dependence V{9) com- 
puted in Fig. [5] for k = 0.12 has not two but four trans- 
porting intervals. The old transporting intervals (^J*^, ^^) 
and {9i , 62) now have non-transporting windows embed- 
ded within them. The borders of these intervals are asso- 
ciated with the depinning of fluxons and they, of course, 
do not exactly coincide with the points of period-doubling 
bifurcations of Fig. [7] The actual depinning takes place 
after the sequence of the period-doubling bifurcations 
and transition to chaotic dynamics. 



C. Dynamical properties of depinned trajectories 

Thus, we have seen that the loss of stability of a pinned 
mode-locked fluxon takes place either due to a tangen- 
tial bifurcation or after a sequence of period-doubling 
bifurcations. In this subsection, we focus on the fluxon 
dynamics at the depinning threshold. 

Consider first the case of mode-locked state stability 
loss via the tangential bifurcation as the phase shift 9 
varies. In this case, we consider the already unstable 
mode-locked state with one of the eigenvalues slightly 
out of the unit circle. After perturbing this state in the 
direction of the unstable eigenvector, the directed fluxon 
propagation begins. It has clear chaotic nature as shown 
in Fig. IsFa). The fluxon dynamics takes place according 




t n 

FIG. 8: Panel (a). Time evolution of central junction phase 
(t>M/2 at K = 0.12, El ^ E2 = 0.27, uj = 0.35, a = 0.1 
for 6 = -0.92462 (1), d = -0.9246 (2), 6 = -0.9245 (3), 
and 9 — —0.924 (4). Panel (b) shows power spectrum for 
9 = -0.924. 

to the intcrmittency type-I scenario: the junction phase 
(in this case the central one with n = N/2) oscillates 
regularly around the equilibrium position as shown by 
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the plateaux of constant phase <j)N/2 Etnd then suddenly 
jumps in a certain direction of propagation (sometimes 
after several chaotic back and forward jumps). The aver- 
age length of the regular (laminar) plateaux is inversely 
proportional to the voltage drop. This length decreases 
as 9 moves away from the bifurcation point. Note the ex- 
tremely sharp dependence of the modulus of the unstable 
eigenvalue |A„„si| as a function of Ei and in Figs.|3]and 
|6j respectively. Since the length of the laminar regime is 
proportional to l/(|A„„st| — 1), it is necessary to change 
the respective parameter in the neighborhood of the de- 
pinning transition with very small increments. Due to 
limited computing power this procedure was not done 
and, as a result, the dependencies in Fig. [5] seem to start 
"out of nowhere" . 

The power spectrum of the junction's phase velocity 
defined as 



m) 



+00 



4>N/2{t)t 



iVtt 



dt 



(9) 



has been computed in Fig. [sjja). The sharp peaks are lo- 
cated at the multiples of the driving frequency: f2 = nw, 
n = 1,2, . . .. These peaks clearly illustrate the remnants 
of the mode-locked dynamics before the tangential bifur- 
cation. The type-I intermittency in the depinning tran- 
sition for the single-harmonically driven kinks has been 
observed in Ref. [4W. 

Now we focus ourselves on the depinning scenario due 
to the period-doubling bifurcations. Consider the case 
depicted in Fig. [7] when the standing fluxon loses its sta- 
bility via the period-doubling bifurcation, for example, 
at K = 0.12, 9 ~ —1.85. The further decrease of 9 leads 
to the sequence of period-doubling bifurcations where 
the initially stable 2T, AT, . . . cycles lose their stabil- 
ity. A typical period-doubling route to chaos takes place 
and the standing fluxon becomes quasiperiodic and later 
chaotic. With decrease of 9, the amplitude of chaotic 
oscillations grows and the fluxon starts to propagate. In 
Fig.joja) the typical trajectories after the depinning point 
are shown. Curve 1 corresponds to non-transporting 
trajectories both regular and chaotic. Since they over- 
lap, the inset with the Poincare sections, computed on 
the plane {4>n/2t4'n/2\ at every time interval f„ — nT, 
n = 1,2, .. . is presented. Curves 2 — 4 correspond to 
transporting trajectories. The dynamical fluxon trajec- 
tories have the structure similar to the case described in 
the previous paragraph. They consist of long intervals, 
where the fluxon stays pinned, which are interrupted by 
short chaotic bursts. During these bursts (which occur 
at random) the fluxon jumps normally in the direction 
specified by the asymmetry of the bias function E{t). 
Sometimes a series of forward and backward jumps takes 
place. Again, the length of the pinned intervals decreases 
as 9 moves away from the critical depinning value. These 
dynamical pictures bear a significant difference, although 
they look similar. In the first case, the standing flux- 
ons are always regular and the transporting trajectory 
is chaotic with long regular intervals where dynamics is 
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FIG. 9: Panel (a). Time evolution of the central junction 
phase 4>M/2 at K = 0.12, Ei = E2 = 0.27, = 0.35 for 
e = -2.08 (1), e = -2.11 (2), e = -2.115 (3) and 9 = -2.12 
(4). The inset shows the Poincare section for 6 — —2.045 (+) 
and 9 = —2.08 (dots). Panels (b)-(d) show Fourier power 
spectra for 9 = -2.05 (b), 9 = -2.055 (c) and 9 = -2.08 (d). 



very close to the periodic one with the period T. In the 
depinning scenario, driven by the period-doubling bifur- 
cations, the chaotization happens before depinning and 
the stationary intervals of the depinned trajectories can- 
not be called laminar because their structure is in fact 
chaotic. 

The power spectra in Fig. [9](b-d) clearly indicate the 
period-doubling route to chaos. Indeed, apart from the 
peaks at nuj, n — 1,2,..., one can clearly see the peaks at 
nuj^mijj /2 and na;±mw/4, to, rt = 1, 2, . . . , in Fig.[9](b). 
As 9 decreases into the chaotic region, the peaks that 
correspond to the aj/4 contribution get smeared out and 
only the uj/2 peaks survive [see panel (c)]. After the 
further decrease of 9 the subharmonic peaks become even 
less pronounced. 



VI. DISCUSSION AND CONCLUSIONS 

In this paper, we have studied the fluxon dynamics in a 
highly discrete annular Josephson junction array driven 
by the asymmetric periodic bias current with the zero 
mean value. It is already well-known that this bias (con- 
sisting of a cosine harmonic and its second overtone) leads 
to a directed fluxon motion which is manifested by the 
non-zero voltage drop. It is interesting to note that for 
the strongly discrete JJA the ratchet transport is mainly 
chaotic, while the non-transporting states are regular. 
Thus, the chaotic dynamics can be identified experimen- 
tally directly from the IV curves. 

We mainly focus on the depinning process of the fluxon 
in the limit of the weak coupling between the neighbor- 
ing junctions (k ^ 1). It appears that the fluxon motion 
is possible in this case provided the amplitude of the ac 
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bias is sufRciently large. The existence diagram on the 
parameter plane {k,Ei) has been computed. On this 
plane, the curve E'^'^'^{k) separates the area where only 
mode-locked standing fluxons exist for any phase shift 6 
from the area where moving fluxons can exist for some 
values of 8. In fact, this diagram can be treated as the 
projection of a more complicated three-dimensional pa- 
rameter space {k,Ei,6) on the plane {k,Ei). 

We have investigated the depinning of the initially 
standing mode-locked fluxon state not far from the crit- 
ical line e['^\k) by analyzing its Floquet spectrum. 
The depinning occurs through chaotization of the mode- 
locked state. The mechanisms of chaotization are diverse 
and depend on the initial position of the non-transporting 
limit cycle in the parameter space {k,Ei,0). The most 
common depinning scenario occurs through the type-I 
intermittency which evolves after the tangential bifurca- 
tion destroys a standing fluxon limit cycle. This scenario 
is ubiquitous when approaching the e['^\k) curve from 
below by increasing the bias amplitude Ei . Another sce- 
nario develops through a sequence of period-doubling bi- 
furcations. It happens slightly above the e['^\k) curve 
when 9 is varied. 

Finally, we would like to point out some further re- 
search directions. The main difficulty in studying the 
kink dynamics in strongly discrete systems is the lack 



of an adequate analytical approximation. As a result, we 
are forced to work with numerical methods. The existing 
approximate theories are based on the perturbations of 
the continuum models where kinks are treated as point 
particles (maybe with the internal mode taken into ac- 
count) in the PN potential. This approach works rela- 
tively well if K > 1 but breaks if k <C 1. An important 
tool that could help to obtain the e['^'' (k) curve analyt- 
ically is the Melnikov criterion. However, the effectively 
use of it requires the explicit expressions for the stable 
and unstable manifolds in order to build the Melnikov 
function. For example, successful implementation of the 
Melnikov method used in Ref. [S31 [53j utilizes the collec- 
tive coordinate approximation which is not applicable in 
our case. Performing the same for the case of small k in 
our opinion is an important challenge. 

In summary, we remark that the biharmonically driven 
ratchet effect in JJAs is robust phenomenon that takes 
place even in the limit of very strong coupling between 
the neighboring junctions. 
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